This notebook is an exercise in the Geospatial Analysis course. You can reference the tutorial at this link.


import math
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon

import folium
from folium import Choropleth, Marker
from folium.plugins import HeatMap, MarkerCluster

from learntools.core import binder
binder.bind(globals())
#from learntools.geospatial.ex5 import *

You'll use the embed_map() function to visualize your maps.

def embed_map(m, file_name):
    from IPython.display import IFrame
    m.save(file_name)
    return IFrame(file_name, width='100%', height='500px')

Exercises

1) Visualize the collision data.

Run the code cell below to load a GeoDataFrame collisions tracking major motor vehicle collisions in 2013-2018.

collisions = gpd.read_file("archive/NYPD_Motor_Vehicle_Collisions/NYPD_Motor_Vehicle_Collisions/NYPD_Motor_Vehicle_Collisions.shp")
collisions.head()
DATE TIME BOROUGH ZIP CODE LATITUDE LONGITUDE LOCATION ON STREET CROSS STRE OFF STREET ... CONTRIBU_2 CONTRIBU_3 CONTRIBU_4 UNIQUE KEY VEHICLE TY VEHICLE _1 VEHICLE _2 VEHICLE _3 VEHICLE _4 geometry
0 07/30/2019 0:00 BRONX 10464 40.841100 -73.784960 (40.8411, -73.78496) None None 121 PILOT STREET ... Unspecified None None 4180045 Sedan Station Wagon/Sport Utility Vehicle Station Wagon/Sport Utility Vehicle None None POINT (1043750.211 245785.815)
1 07/30/2019 0:10 QUEENS 11423 40.710827 -73.770660 (40.710827, -73.77066) JAMAICA AVENUE 188 STREET None ... None None None 4180007 Sedan Sedan None None None POINT (1047831.185 198333.171)
2 07/30/2019 0:25 None None 40.880318 -73.841286 (40.880318, -73.841286) BOSTON ROAD None None ... None None None 4179575 Sedan Station Wagon/Sport Utility Vehicle None None None POINT (1028139.293 260041.178)
3 07/30/2019 0:35 MANHATTAN 10036 40.756744 -73.984590 (40.756744, -73.98459) None None 155 WEST 44 STREET ... None None None 4179544 Box Truck Station Wagon/Sport Utility Vehicle None None None POINT (988519.261 214979.320)
4 07/30/2019 10:00 BROOKLYN 11223 40.600090 -73.965910 (40.60009, -73.96591) AVENUE T OCEAN PARKWAY None ... None None None 4180660 Station Wagon/Sport Utility Vehicle Bike None None None POINT (993716.669 157907.212)

5 rows × 30 columns

Use the "LATITUDE" and "LONGITUDE" columns to create an interactive map to visualize the collision data. What type of map do you think is most effective?

m_1 = folium.Map(location=[40.7, -74], zoom_start=11) 

# Your code here: Visualize the collision data
HeatMap(data=collisions[['LATITUDE', 'LONGITUDE']], radius=15).add_to(m_1)

# Uncomment to see a hint
#q_1.hint()
m_1
# Show the map
#embed_map(m_1, "q_1.html")
Make this Notebook Trusted to load map: File -> Trust Notebook
q_1.check()

# Uncomment to see our solution (your code may look different!)
#q_1.solution()

2) Understand hospital coverage.

Run the next code cell to load the hospital data.

hospitals = gpd.read_file("archive/nyu_2451_34494/nyu_2451_34494/nyu_2451_34494.shp")
hospitals.head()
id name address zip factype facname capacity capname bcode xcoord ycoord latitude longitude geometry
0 317000001H1178 BRONX-LEBANON HOSPITAL CENTER - CONCOURSE DIVI... 1650 Grand Concourse 10457 3102 Hospital 415 Beds 36005 1008872.0 246596.0 40.843490 -73.911010 POINT (1008872.000 246596.000)
1 317000001H1164 BRONX-LEBANON HOSPITAL CENTER - FULTON DIVISION 1276 Fulton Ave 10456 3102 Hospital 164 Beds 36005 1011044.0 242204.0 40.831429 -73.903178 POINT (1011044.000 242204.000)
2 317000011H1175 CALVARY HOSPITAL INC 1740-70 Eastchester Rd 10461 3102 Hospital 225 Beds 36005 1027505.0 248287.0 40.848060 -73.843656 POINT (1027505.000 248287.000)
3 317000002H1165 JACOBI MEDICAL CENTER 1400 Pelham Pkwy 10461 3102 Hospital 457 Beds 36005 1027042.0 251065.0 40.855687 -73.845311 POINT (1027042.000 251065.000)
4 317000008H1172 LINCOLN MEDICAL & MENTAL HEALTH CENTER 234 E 149 St 10451 3102 Hospital 362 Beds 36005 1005154.0 236853.0 40.816758 -73.924478 POINT (1005154.000 236853.000)

Use the "latitude" and "longitude" columns to visualize the hospital locations.

m_2 = folium.Map(location=[40.7, -74], zoom_start=11) 

# Your code here: Visualize the hospital locations
for idx, row in hospitals.iterrows():
    Marker([row['latitude'], row['longitude']]).add_to(m_2)

# Uncomment to see a hint
#q_2.hint()
m_2
# Show the map
#embed_map(m_2, "q_2.html")
Make this Notebook Trusted to load map: File -> Trust Notebook
q_2.check()

# Uncomment to see our solution (your code may look different!)
#q_2.solution()

3) When was the closest hospital more than 10 kilometers away?

Create a DataFrame outside_range containing all rows from collisions with crashes that occurred more than 10 kilometers from the closest hospital.

Note that both hospitals and collisions have EPSG 2263 as the coordinate reference system, and EPSG 2263 has units of meters.

ten_kilom_buffer = hospitals.geometry.buffer(10*1000)
ten_kilom_buffer.head()
union=ten_kilom_buffer.geometry.unary_union
li=[]
for idx, row in collisions.iterrows():
    if(union.contains(row.geometry)) is False:
                                         li.append(idx)
outside_range = collisions.iloc[li]
outside_range.head()

# Check your answer
#q_3.check()
DATE TIME BOROUGH ZIP CODE LATITUDE LONGITUDE LOCATION ON STREET CROSS STRE OFF STREET ... CONTRIBU_2 CONTRIBU_3 CONTRIBU_4 UNIQUE KEY VEHICLE TY VEHICLE _1 VEHICLE _2 VEHICLE _3 VEHICLE _4 geometry
0 07/30/2019 0:00 BRONX 10464 40.841100 -73.784960 (40.8411, -73.78496) None None 121 PILOT STREET ... Unspecified None None 4180045 Sedan Station Wagon/Sport Utility Vehicle Station Wagon/Sport Utility Vehicle None None POINT (1043750.211 245785.815)
1 07/30/2019 0:10 QUEENS 11423 40.710827 -73.770660 (40.710827, -73.77066) JAMAICA AVENUE 188 STREET None ... None None None 4180007 Sedan Sedan None None None POINT (1047831.185 198333.171)
5 07/30/2019 10:50 QUEENS 11423 40.721060 -73.759450 (40.72106, -73.75945) FRANCIS LEWIS BOULEVARD HILLSIDE AVENUE None ... None None None 4179812 Sedan Box Truck None None None POINT (1050928.749 202069.687)
6 07/30/2019 10:55 QUEENS 11434 40.676228 -73.761120 (40.676228, -73.76112) CRANDALL AVENUE CHENEY STREET None ... None None None 4180464 Station Wagon/Sport Utility Vehicle Station Wagon/Sport Utility Vehicle None None None POINT (1050510.380 185734.852)
15 07/30/2019 13:05 None None 40.588413 -74.166725 (40.588413, -74.166725) None None 26 RICHMOND HILL ROAD ... None None None 4180091 Station Wagon/Sport Utility Vehicle Box Truck None None None POINT (937943.004 153695.210)

5 rows × 30 columns

#q_3.hint()
#q_3.solution()
len(li)
39595

The next code cell calculates the percentage of collisions that occurred more than 10 kilometers away from the closest hospital.

percentage = round(100*len(outside_range)/len(collisions), 2)
print("Percentage of collisions more than 10 km away from the closest hospital: {}%".format(percentage))
Percentage of collisions more than 10 km away from the closest hospital: 15.12%

4) Make a recommender.

When collisions occur in distant locations, it becomes even more vital that injured persons are transported to the nearest available hospital.

With this in mind, you decide to create a recommender that:

  • takes the location of the crash (in EPSG 2263) as input,
  • finds the closest hospital (where distance calculations are done in EPSG 2263), and
  • returns the name of the closest hospital.
def best_hospital(collision_location):
    # Your code here
    distances = hospitals.geometry.distance(collision_location)

    name = hospitals.iloc[distances.idxmin()]['name']
    return name

# Test your function: this should suggest CALVARY HOSPITAL INC
print(best_hospital(outside_range.geometry.iloc[0]))

# Check your answer
#q_4.check()
CALVARY HOSPITAL INC
#q_4.hint()
#q_4.solution()

5) Which hospital is under the highest demand?

Considering only collisions in the outside_range DataFrame, which hospital is most recommended?

Your answer should be a Python string that exactly matches the name of the hospital returned by the function you created in 4).

highest_demand = "JAMAICA HOSPITAL MEDICAL CENTER"

# Check your answer
#q_5.check()
#q_5.hint()
#q_5.solution()

6) Where should the city construct new hospitals?

Run the next code cell (without changes) to visualize hospital locations, in addition to collisions that occurred more than 10 kilometers away from the closest hospital.

m_6 = folium.Map(location=[40.7, -74], zoom_start=11) 


folium.GeoJson(ten_kilom_buffer.geometry.to_crs(epsg=4326)).add_to(m_6)
HeatMap(data=outside_range[['LATITUDE', 'LONGITUDE']], radius=9).add_to(m_6)
folium.LatLngPopup().add_to(m_6)
for idx, rows in hospitals.iterrows():
    Marker([rows['latitude'], rows['longitude']], popup=rows['name']).add_to(m_6)
folium.LatLngPopup().add_to(m_6)


m_6
#embed_map(m_6, 'm_6.html')
Make this Notebook Trusted to load map: File -> Trust Notebook

Click anywhere on the map to see a pop-up with the corresponding location in latitude and longitude.

The city of New York reaches out to you for help with deciding locations for two brand new hospitals. They specifically want your help with identifying locations to bring the calculated percentage from step 3) to less than ten percent. Using the map (and without worrying about zoning laws or what potential buildings would have to be removed in order to build the hospitals), can you identify two locations that would help the city accomplish this goal?

Put the proposed latitude and longitude for hospital 1 in lat_1 and long_1, respectively. (Likewise for hospital 2.)

Then, run the rest of the cell as-is to see the effect of the new hospitals. Your answer will be marked correct, if the two new hospitals bring the percentage to less than ten percent.

lat_1 = 40.6728
long_1 = -73.7478

# Your answer here: proposed location of hospital 2
lat_2 = 40.6791
long_2 = -73.8666


# Do not modify the code below this line
try:
    new_df = pd.DataFrame(
        {'Latitude': [lat_1, lat_2],
         'Longitude': [long_1, long_2]})
    new_gdf = gpd.GeoDataFrame(new_df, geometry=gpd.points_from_xy(new_df.Longitude, new_df.Latitude))
    new_gdf.crs = {'init' :'epsg:4326'}
    new_gdf = new_gdf.to_crs(epsg=2263)
    # get new percentage
    new_coverage = gpd.GeoDataFrame(geometry=new_gdf.geometry).buffer(10000)
    new_my_union = new_coverage.geometry.unary_union
    new_outside_range = outside_range.loc[~outside_range["geometry"].apply(lambda x: new_my_union.contains(x))]
    new_percentage = round(100*len(new_outside_range)/len(collisions), 2)
    print("(NEW) Percentage of collisions more than 10 km away from the closest hospital: {}%".format(new_percentage))
    
    #q_6.check()
    
    # make the map
    m = folium.Map(location=[40.7, -74], zoom_start=11) 
    folium.GeoJson(coverage.geometry.to_crs(epsg=4326)).add_to(m)
    folium.GeoJson(new_coverage.geometry.to_crs(epsg=4326)).add_to(m)
    for idx, row in new_gdf.iterrows():
        Marker([row['Latitude'], row['Longitude']]).add_to(m)
    HeatMap(data=new_outside_range[['LATITUDE', 'LONGITUDE']], radius=9).add_to(m)
    folium.LatLngPopup().add_to(m)
    display(embed_map(m, 'q_6.html'))
except:
    q_6.hint()
c:\Users\lg\anaconda3\lib\site-packages\pyproj\crs\crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  in_crs_string = _prepare_from_proj_string(in_crs_string)
(NEW) Percentage of collisions more than 10 km away from the closest hospital: 9.12%
#q_6.solution()

Congratulations!

You have just completed the Geospatial Analysis micro-course! Great job!


Have questions or comments? Visit the course discussion forum to chat with other learners.